Tarea 4: Bayesian Inference Part I

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Fecha límite de entrega:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Primera Parte: Preguntas Teóricas
  5. Segunda Parte: Elaboración de Código

Objetivo

a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

Primera Parte: Preguntas Teóricas

A continuación, se presentaran diferentes preguntas que abordan las temáticas vistas en clases. Por favor responda cada una de estas preguntas de forma breve, no más de 4 o 5 lineas.

Pregunta 1:

Explique cual es la diferencia fundamental entre la estadística bayesiana y frecuentista.

La diferencia fundamental entre la estadística Bayesiana con la Frecuentista es la misma definición de probabilidad. En el mundo Frecuentista estamos acostumbrados a pensar la probabilidad de un suceso como un valor definido, al cuál se tiende cuando estamos en medidas a gran escala (la idea base de la Ley de los Grandes Números) y el cuál evita toda subjetividad, pero en el mundo Bayesiano la probabilidad es una creencia o convicción de que ocurra un suceso, y este depende de la información a priori que se tiene del evento.

Pregunta 2:

Discuta la siguiente afirmación La inferencia bayesiana permite fácilmente utilizar distintos tipos de información.

Debido a que la inferencia Bayesiana necesita un prior para poder generar una distribución del posterior, es posible adaptar esta información inicial basado en información externa que permita entender de mejor forma el comportamiento del parámetro. En otras palabras, es posible utilizar información externa respecto al posible comportamiento del parámetro, algo imposible en la inferencia frecuentista que trabajaba directamente con los datos.

Pregunta 3:

Explique la diferencia entre prior probability y posterior probability

La Prior Probability se refieren a las creencias iniciales respecto a un parámetro, mientras que la Posterior Probability son las creencias actualizadas utilizando los datos y lo que utilizamos para generar una distribución. A medida que recibimos mayor cantidad de datos, los que fueron en algún momento Posterior pasan a ser la creencia inicial con el fin de actualizar el modelo.

Pregunta 4:

El estadista bayesiano “Bruno Finetti” menciona la siguiente frase en su libro de probabilidad: La probabilidad no existe. Lo que en verdad quizo decir es que la probabilidad es un método para describir incertidumbre en un observador con conocimiento limitado. Discuta esta información utilizando el ejemplo del lanzamiento del globo terraqueo visto en clases. ¿Que significa decir “la probabilidad de que sea agua es 0.7”?

La frase mencionada funciona bien y se entiende en el mundo frecuentista, ya que se refiere a que un parámetro desconocido en un principio tiene un valor exacto y real, que en este caso es la proporción Agua-Tierra en el planeta. Pero esto no es así en el mundo Bayesiano, ya que el parámetro es una variable aleatoria por sí sola, por lo cuál este contiene cierto margen de error e incertidumbre imposibles de reducir a 0.

Pregunta 5:

¿ Que ventaja entrega que la distribución de la posterior este en la misma familia de distribución de probabilidad que la del prior?. De un ejemplo de alguna distribución que posee este comportamiento.

La ventaja principal se debe a la utilización del método conocido como Priors Conjugados. En resumen, conocer la distribución del prior y del posterior ayuda a modelar el problema de forma analítica y a su vez de forma mucho menos costosa computacionalmente, debido a que en vez de calcular \(f(d)\) y la Likelihood por separado para pasar a computar el posterior, podemos tomar un atajo y resolver directamente el posterior.

Un ejemplo es cuando las distribuciones siguen una distribución Beta. En este caso, en vez de resolver el prior (Beta) y likelihood para calcular el posterior, podemos calcular directamente el posterior como una función Beta con valores distintos.

Pregunta 6:

Señale y explique los pasos de la grid approximation para obtener la posterior y responda las siguientes preguntas:

  1. ¿Cual de los pasos señalados nos permite obtener una distribución de la posterior mas precisa?.
  2. Cuales son las limitaciones de la grid approximation.

Primero, se define la grilla con la cantidad de puntos que queremos utilizar. A mazor cantidad de puntos, más exacta y precisa es la distribución (parte a),

Segundo, se computa el valor del prior por cada valor posible de la grilla (\(f(\theta)\))

Tercero, Se computa el Likelihood por cada valor posible del parámetro, (\(f(d|\theta)\))

Cuarto, Se computa el posterior sin estandarizar por cada valor posible del parámetro (\(f(d|\theta) f(\theta)\))

Finalmente estandarizamos el posterior al dividir por la suma de los valores (\(/ f(d)\)).

La principal limitación de la grid approximation es que es poco escalable, ya que a grandes cantidades de datos este es más complejo computacionalmente debido a la cantidad de operaciones a realizar.

Pregunta 7:

¿ Por qué es necesario aprender a trabajar con muestras de la posterior?.

Porque así podemos modelar correctamente la distribución de valores posibles que hayan generado los datos, y ajustar el modelo a lo esperado. Por ejemplo, si no trabajamos correctamente con las muestras de la posterior, podríamos obtener una distribución que no contenga el valor real como parte de los valores más probables, lo que invalida la credibilidad y calidad del modelo.

Pregunta 8:

Señale si las siguientes afirmaciones son verdaderas o falsas, en el caso que sean falsas justifique su respuesta:

  • Los point estimates de la posterior no entregan información relevante en un estudio.

  • Un intervalo de confianza es un intervalo dentro del cual un valor de parámetro no observado cae con una probabilidad particular, mientras que un intervalo de credibilidad es un rango de valores en el que se estima que estará cierto valor desconocido.

  • La principal ventaja de HPDI frente a un intervalo de credibilidad es que si la posterior no distribuye de forma normal, el HPDI será capaz de detectar los puntos de interés, mientras que un intervalo de credibilidad lo ignoraría al asumir simetría.

Falso. Cuando la posterior distribuye de forma normal, los point estimates son una forma alternativa de resumir los valores más relevantes como la moda o mediana. Sin embargo, estos no son tan útiles como los intervalos de confianza, pero no significa que no entreguen información relevante.

Falso. La descripción indicada corresponde a un intervalo de confianza en el mundo Frecuentista, no Bayesiano. Los intervalos de credibilidad no se refieren a la probabilidad de contener el valor real, si no que representa la creencia de valores que puede tomar el estadístico muestral.

Verdadero.

Pregunta 9

Suponga que tiene dos especies de pandas. Cada una de las especies es igual de común y es imposible distinguirlas físicamente. Una de las diferencias entre las especies es el tamaño de sus familias. Si denotamos por \(\theta\) a la especie del panda se tiene que, cuando la especie es \(\theta = 1\) tiene dos bebes un \(10\%\) de las veces mientras que la especie \(\theta = 2\) tiene dos bebes un \(20\%\) de las veces, mientras que el resto de veces ambas especies tienen un solo bebe.

Suponga que usted esta intentando determinar la especie de un panda que que tiene como registro de nacimientos al conjunto \(D\), considere quela especie de un panda que acaba de dar a luz a dos bebes, es decir \(D = \{\text{dos bebes}\}\)

  • ¿Cual es la probabilidad de que pertenezca a la especie \(1\)?

  • Suponga ahora que el mismo panda acaba de dar a luz y esta vez es solo un bebe. Calcule la probabilidad posterior de que el panda sea de especie \(1\). ¿Que cambia con el calculo anterior?

  • Suponga que le ofrecen hacer un test genético a su panda, como suele ser común con los test no es perfecto y le mencionan las siguientes características:

    • La probabilidad de que correctamente idenfitique a la especie \(1\) es de \(0.8\)
    • La probabilidad de que correctamente identifique a la especie \(2\) es de \(0.65\)

Se administra el test y se obtiene un resultado positivo a la especie \(1\). Sin utilizar la información en \(D\) calcule la probabilidad posterior de que su panda sea de la especie \(1\). Repita sus cálculos utilizando la información recopilada en \(D\). ¿En que varían sus resultados?

Idk

Segunda Parte: Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.0.5
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.1.0     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.0.5
## Warning: package 'tibble' was built under R version 4.0.5
## Warning: package 'tidyr' was built under R version 4.0.5
## Warning: package 'readr' was built under R version 4.0.5
## Warning: package 'dplyr' was built under R version 4.0.5
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(dplyr)
library(tidyr)
library(purrr)

# Para realizar plots
#library(scatterplot3d)
library(ggplot2)
library(plotly)
## Warning: package 'plotly' was built under R version 4.0.5
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.0.5
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(grid)

# Análisis bayesiano
library(rethinking)
## Loading required package: rstan
## Warning: package 'rstan' was built under R version 4.0.5
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 4.0.5
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:tidyr':
## 
##     extract
## Loading required package: parallel
## rethinking (Version 2.13)
## 
## Attaching package: 'rethinking'
## The following object is masked from 'package:purrr':
## 
##     map
## The following object is masked from 'package:stats':
## 
##     rstudent

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Preguntas: Introducción a Grid Approximation

Las primeras dos preguntas de esta tarea tienen como objetivo introducirlos en la inferencia Bayesiana utilizando la técnica Grid Approximation para obtener una aproximación de la posterior. Al finalizar los problemas ustedes deberán ser capaces de visualizar los efectos que tiene el prior en la posterior, saber cómo realizar una Grid Approximation y comprender como utilizar Percentile Interval (PI) en una posterior.

Pregunta 1:

Considere el dataset “moneda.csv” donde se encuentran los resultados de un experimento lanzando una moneda, el objetivo de esta pregunta es estudiar mediante técnicas de inferencia Bayesiana el valor de la probabilidad de que salga cara, representado por el valor \(1\). Puede usar la librerira rethinking durante toda esta pregunta (si lo desea).

dataMoneda <- read.csv("moneda.csv", header = TRUE)

  • Programe el metodo grid approximation para distintos tamaños de experimento. ¿Como van cambiando las curvas posterior?

  • Repita el mismo análisis anterior pero utilizando el método de Laplace (no necesita programar el método, puede utilizar la libreria “rethinking”). ¿Como se comparan con los resultados anteriores?.

  • Grafique la densidad de la posterior y encuentre la proporción de los siguientes defined boundaries:

    • \((0, 0.4)\)
    • \((0.4,0.7)\)
    • \((0.7,1)\)

¿Como puede interpretar los resultados?

  • Calcule un intervalo de credibilidad al \(50\%\), \(75\%\) y \(95\%\). ¿Como se puede interpretar los resultados? ¿Cual podría ser un problema al usar intervalos de credibilidad?.
  • Genere los intervalos HDPI para \(50\%\), \(75\%\) y \(95\%\), compárelos con los intervalos de credibilidad de la parte anterior. ¿En que se diferencian los intervalos de credibilidad con los HDPI?.

Respuesta

Aproximación de Grilla

# Función de grid approximation, que recibe un data frame y la cantidad de valores en la grilla
grid_app <- function(data, length){
  # Generamos la grilla del tamaño requerido y con un prior cte
  p_grid <- seq(from=0 , to=1 , length.out=length)
  prior <- rep(1 , length)
  
  # Contamos los éxitos registrados en el dataframe
  count <- sum(data[data$Resultado == 1,])
  # Generamos la likelihood como una distribución binomial con la cantidad de exitos registrada
  likelihood <- dbinom(count , size=1000 , prob=p_grid)
  
  # Finalmente calculamos la posterior y se plotea
  posterior <- likelihood * prior
  posterior <- posterior / sum(posterior)
  
  title <- paste("size =", length)
  
  data <- data.frame(posterior = posterior, grid = p_grid)
  row <- data[which.max(data$posterior),]
  
  #Printeamos el valor máximo de la posterior
  print(row)
  
  plot <- ggplot(data, aes(x=grid, y=posterior)) +
    geom_line() + 
    geom_vline(xintercept = row$grid, color = "red", size=0.5) +
    xlab("grid") +  
    ylab("posterior") + 
    ggtitle(title) + 
    theme(plot.title = element_text(hjust = 0.5))
  
  # Graficamos además la densidad
  samples <- sample(p_grid , prob=posterior , size=1e4 , replace=TRUE)
  dens(samples)
  
  return(list(plot=plot, sample = samples))
}
plot1 <- grid_app(dataMoneda, 5)[[1]]
##   posterior grid
## 3         1  0.5

plot2 <- grid_app(dataMoneda, 10)[[1]]
##   posterior      grid
## 6         1 0.5555556

plot3 <- grid_app(dataMoneda, 100)[[1]]
##    posterior      grid
## 56 0.2499644 0.5555556

plot4 <- grid_app(dataMoneda, 1000)[[1]]
##      posterior      grid
## 552 0.02540291 0.5515516

grid.arrange(plot1, plot2, plot3, plot4,
             # Numero de filas y columnas en los que los queremos sortear
             ncol=2, nrow=2,
             # Títulos
             top = textGrob("Grid Approximation of dataMoneda, using different grid sizes", 
                            gp=gpar(fontsize=15)))

Es claro que a medida que aumentamos la cantidad de puntos en la grilla, tenemos una visualización más exacta de la posterior. La última gráfica muestra una comparación de los 4 experimentos con el tamaño de la grilla creciente, donde se observa que cada vez se obtiene una distribución con menor desviación estándar y con una mediana más exacta. Estas visualizaciones están más detalladas en los print de los samples de cada experimento.

Aproximación de Laplace

laplace_app <- function(data){
  # Contamos los aciertos
  count <- sum(data[data$Resultado == 1,])
  
  # Utilizamos la funcion quap del paquete rethinking
  coin.qa <- quap(
  alist(
    C ~ dbinom(C+S, p) , # binomial likelihood
    p ~ dunif(0,1)   # uniform prior
  ) ,
  data=list(C=count, S=(1000-count)))

  # display summary of quadratic approximation
  print(precis(coin.qa))
  sample.quap <- extract.samples(coin.qa)
  return(sample.quap)
}
samples <- laplace_app(dataMoneda)
##    mean         sd      5.5%     94.5%
## p 0.552 0.01572558 0.5268675 0.5771325
dens(samples)

Notemos que la mediana sigue estando al rededor del valor 0.55, y que la forma de la distribución de la poterior es muy similar al gráfico de densidad de la aproximación con grilla de tamaño 1000. Esto quiere decir que para este caso, ambas aproximaciones se comportan de buena forma y arrojan resultados muy similares

Gráfico de densidad y proporciones:

Grafiquemos primero la densidad

# Utilizamos los samples retornados por la aproximacion de grilla
samples <- grid_app(dataMoneda, 1000)[[2]]
##      posterior      grid
## 552 0.02540291 0.5515516

# Generamos un dataframe con los datos
laplace <- data.frame(p=samples)
p1 <- ggplot(laplace, aes(x=p)) + 
  geom_density(color="black", fill = "red", alpha=0.3) + 
  # Graficamos la densidad y la mediana registrada
  geom_vline(xintercept = 0.5519999, color = "red", size=1) +
  xlim(0.5, 0.61) +
  xlab("p") +  
  ylab("posterior") + 
  ggtitle("Density plot") + 
  theme(plot.title = element_text(hjust = 0.5))
p1
## Warning: Removed 7 rows containing non-finite values (stat_density).

# Calculamos ahora la proporcion por cada sector solicitao
prop1 = sum(laplace[laplace$p < 0.4, ])
prop2 = sum(laplace[laplace$p > 0.4 && laplace$p < 0.7, ])
prop3 = sum(laplace[laplace$p > 0.7, ])

prop1/sum(laplace$p)
## [1] 0
prop2/sum(laplace$p)
## [1] 1
prop3/sum(laplace$p)
## [1] 0

Debido a que ambos lados contienen el mismo porcentaje de la densidad del posterior, tenemos que estamos frente a una distribución simétrica, y totalmente centrada en los rangos del 0,4 al 0,7.

Intervalos de credibilidad

Primer intervalo, 50%:

# Utilizamos la funcion cuantia para calcular los limites 
prop1 = quantile(samples, c(0.25, 0.75))
prop1
##       25%       75% 
## 0.5415415 0.5625626
# Generamos un primer plot que contenga la densidad, los dos limites
# calculados de la parte anterior, y la mediana de la posterior
p1 <- ggplot(laplace, aes(x=p)) + 
  geom_density(color="black") + 
  geom_vline(xintercept = 0.5415415, color = "red", size=1) +
  geom_vline(xintercept = 0.5625626, color = "red", size=1) + 
  geom_vline(xintercept = 0.5519999, color = "blue", size=1)
  
# Generamos un build para hacer fill del area encontrada
d = ggplot_build(p1)$data[[1]]

# Llenamos el area solicitada y añadimos los titulos y limites
p1 <- p1 + 
  # Codigo que hace fill de la area calculada
  geom_area(data = subset(d, x>0.5415415 & x<0.5625626), aes(x=x,y=y), fill = "red", alpha = 0.5) +
  xlim(0.5, 0.61) +
  xlab("p") +  
  ylab("posterior") + 
  ggtitle("Density plot, 50% centered") + 
  theme(plot.title = element_text(hjust = 0.5))
p1
## Warning: Removed 7 rows containing non-finite values (stat_density).

Segundo intervalo, 75%:

prop2 = quantile(samples, c(0.125, 0.875))
prop2
##     12.5%     87.5% 
## 0.5335335 0.5705706
p2 <- ggplot(laplace, aes(x=p)) + 
  geom_density(color="black") + 
  geom_vline(xintercept = 0.5335335, color = "red", size=1) +
  geom_vline(xintercept = 0.5705706, color = "red", size=1) +  
  geom_vline(xintercept = 0.5519999, color = "blue", size=1)
  
d = ggplot_build(p2)$data[[1]]

p2 <- p2 + 
  geom_area(data = subset(d, x>0.5335335 & x<0.5705706), aes(x=x,y=y), fill = "red", alpha = 0.5) +
  xlim(0.5, 0.61) +
  xlab("p") +  
  ylab("posterior") + 
  ggtitle("Density plot, 75% centered") + 
  theme(plot.title = element_text(hjust = 0.5))
p2
## Warning: Removed 7 rows containing non-finite values (stat_density).

Tercer intervalo, 95%:

prop3 = quantile(samples, c(0.025, 0.975))
prop3
##      2.5%     97.5% 
## 0.5215215 0.5835836
p3 <- ggplot(laplace, aes(x=p)) + 
  geom_density(color="black") + 
  geom_vline(xintercept = 0.5205205, color = "red", size=1) +
  geom_vline(xintercept = 0.5825826, color = "red", size=1) +  
  geom_vline(xintercept = 0.5519999, color = "blue", size=1)
  
d = ggplot_build(p2)$data[[1]]
## Warning: Removed 7 rows containing non-finite values (stat_density).
p3 <- p3 + 
  geom_area(data = subset(d, x>0.5205205 & x<0.5825826), aes(x=x,y=y), fill = "red", alpha = 0.5) +
  xlim(0.5, 0.61) +
  xlab("p") +  
  ylab("posterior") + 
  ggtitle("Density plot, 95% centered") + 
  theme(plot.title = element_text(hjust = 0.5))
p3
## Warning: Removed 7 rows containing non-finite values (stat_density).

En este caso la mediana se encuentra siempre en el centro de los límites calculados. Esto sucede debido al comportamiento binomial que tienen las monedas, y a la forma simétrica que presenta la posterior.

Intervalos de credibilidad con HPDI

Primer intervalo, 50%:

prop1HPDI = HPDI(samples, prob=0.5)
prop1HPDI
##      |0.5      0.5| 
## 0.5405405 0.5615616
p1HPDI <- ggplot(laplace, aes(x=p)) + 
  geom_density(color="black") + 
  geom_vline(xintercept = 0.5405405, color = "red", size=1) +
  geom_vline(xintercept = 0.5615616, color = "red", size=1) +  
  geom_vline(xintercept = 0.5519999, color = "blue", size=1)
  
d = ggplot_build(p1HPDI)$data[[1]]

p1HPDI <- p1HPDI + 
  geom_area(data = subset(d, x>0.5405405 & x<0.5615616), aes(x=x,y=y), fill = "red", alpha = 0.5) +
  xlim(0.5, 0.61) +
  xlab("p") +  
  ylab("posterior") + 
  ggtitle("Density plot, 50% HPDI") + 
  theme(plot.title = element_text(hjust = 0.5))
p1HPDI
## Warning: Removed 7 rows containing non-finite values (stat_density).

Segundo intervalo, 75%:

prop2HPDI = HPDI(samples, prob=0.75)
prop2HPDI
##     |0.75     0.75| 
## 0.5325325 0.5685686
p2HPDI <- ggplot(laplace, aes(x=p)) + 
  geom_density(color="black") + 
  geom_vline(xintercept = 0.5315315, color = "red", size=1) +
  geom_vline(xintercept = 0.5675676, color = "red", size=1) +  
  geom_vline(xintercept = 0.5519999, color = "blue", size=1)
  
d = ggplot_build(p2HPDI)$data[[1]]

p2HPDI <- p2HPDI + 
  geom_area(data = subset(d, x>0.5315315 & x<0.5675676), aes(x=x,y=y), fill = "red", alpha = 0.5) +
  xlim(0.5, 0.61) +
  xlab("p") +  
  ylab("posterior") + 
  ggtitle("Density plot, 75% HPDI") + 
  theme(plot.title = element_text(hjust = 0.5))
p2HPDI
## Warning: Removed 7 rows containing non-finite values (stat_density).

Tercer intervalo, 95%:

prop3HPDI = HPDI(samples, prob=0.95)
prop3HPDI
##     |0.95     0.95| 
## 0.5205205 0.5825826
p3HPDI <- ggplot(laplace, aes(x=p)) + 
  geom_density(color="black") + 
  geom_vline(xintercept = 0.5205205, color = "red", size=1) +
  geom_vline(xintercept = 0.5815816, color = "red", size=1) +  
  geom_vline(xintercept = 0.5519999, color = "blue", size=1)
  
d = ggplot_build(p2HPDI)$data[[1]]
## Warning: Removed 7 rows containing non-finite values (stat_density).
p3HPDI <- p3HPDI + 
  geom_area(data = subset(d, x>0.5205205 & x<0.5815816), aes(x=x,y=y), fill = "red", alpha = 0.5) +
  xlim(0.5, 0.61) +
  xlab("p") +  
  ylab("posterior") + 
  ggtitle("Density plot, 95% HPDI") + 
  theme(plot.title = element_text(hjust = 0.5))
p3HPDI
## Warning: Removed 7 rows containing non-finite values (stat_density).

Esta forma alternativa de encontrar los intervalos de confianza no solo obliga que el área encontrada contenga a la mediana (lo cuál debido a la simetría del problema, los intervalos sin HPDI también la contienen), si no que encuentra el área con los limites más cercanos entre sí. Esto hace que ahora la mediana no se encuentre necesariamente al centro del área o intervalo de confianza, y esto puede observarse en los primeros gráficos (en especial el de 75%), donde la mediana se encuentra más a la derecha del centro simétrico del área.

Lo anterior indica que si bien la distribución es simétrica, las curvaturas hacia ambos extremos de la binomial no lo son. Otra forma de verlo es que la distribución está más concentrada hacia la izquierda del valor máximo, y que valores mayores a la media son menos probables que los menores.

Pregunta 2: Grid Approximation y Posterior Predictive Distribution

El objetivo de esta pregunta es comprender el concepto de sample prediction visto en clases y realizar predicciones en base a una posterior.

Un conjunto de carteros aburridos de las mordidas de perros ha decidido realizar un catastro de mordidas recibidas por los empleados de su empresa en un periodo de dos meses, planeando en base a estos datos realizar inferencia bayesiana. Los datos de las mordidas estas datos por el dataset no+mordidas.csv, en donde cada fila representa las mordidas recibidas por diferentes carteros y las columnas señalan si fue mordido o no el cartero en los meses de estudio (notar que si fue mordido sera señalado con un 1, de lo contrario es señalado con un 0). Cabe señalar que un cartero no puede ser mordido mas de una vez al mes, ya que el damnificado recibe licencia por todos los días restantes del mes tras la mordida, reincorporándose el siguiente mes al trabajo.

df = read.csv("no+mordidas.csv")
head(df)

En base a los datos, realice los siguientes puntos:

  • Realice una grid approximation para estimar la probabilidad que un cartero sea mordido, para esto junte los datos del mes 1 y 2 de estudio. Señale el máximo valor de la posterior.

  • Utilizando la posterior obtenida en el paso anterior, utilice rbinom para simular 10.000 réplicas de 500 registros de mordidas. Con esto, deberá obtener 10.000 números, cada uno de los cuales es un recuento de las mordidas obtenidas en el registro de datos. Compare la distribución del número de los carteros mordidos predichos con el número real de los datos (248 carteros mordidos de un total de 500 datos). ¿El modelo se ajusta bien a los datos? Es decir, ¿la distribución de las predicciones incluye la observación real como resultado central y probable?

  • Como se comento en el comienzo bites_month1 contiene las mordidas señaladas por un conjunto de personas en el primer mes. Haciendo uso de bites_month2, obtenga la posterior de que una persona que fue mordida en el primer mes, sea mordida nuevamente en el segundo mes. Para esto, se recomienda comenzar buscando los carteros que fueron mordidos el primer mes y en base a estos generar una búsqueda indexada para obtener el número solicitado. Hecho esto, simule ese número carteros mordidos 10.000 veces. De los resultados obtenidos, compare el recuento de carteros mordidos con el recuento real. ¿Cómo se ve el modelo desde este punto de vista?

Respuesta Aqui

Primera parte, Grid approximation

# Inicializamos la cantidad de valores de la grilla
length <- 1000

# Creamos la grilla con probabilidades iguales
p_grid <- seq(from=0 , to=1 , length.out=length)
prior <- rep(1 , length)

# Contamos los éxitos registrados en el dataframe
count <- sum(df$bites_month_1) + sum(df$bites_month_2)

# Generamos la likelihood como una distribución binomial con la cantidad de exitos registrada
likelihood <- dbinom(count , size=500 , prob=p_grid)

# Finalmente calculamos la posterior y se plotea
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)

# Caculo del maximo valor
data <- data.frame(posterior = posterior, grid = p_grid)
row <- data[which.max(data$posterior),]

plot <- ggplot(data, aes(x=grid, y=posterior)) +
  geom_line() + 
  geom_vline(xintercept = row$grid, color = "red", size=0.5) +
  xlab("grid") +  
  ylab("posterior") + 
  ggtitle("Grid Approximation to no+mordidas dataset") + 
  theme(plot.title = element_text(hjust = 0.5))
plot

# Graficamos además la densidad
samples <- sample(p_grid , prob=posterior , size=1e4 , replace=TRUE)

#Printeamos el valor máximo de la posterior, la media y la densidad
print(row)
##      posterior      grid
## 497 0.01788204 0.4964965
mean(samples)
## [1] 0.4963612
dens(samples)

Obteniendo que el maximo es 0.4964965 y a mean es 0.4962622

Segunda parte, Posterior

# Posterior Predictive
new <- rbinom(1e4, size=500, prob=samples)
dens(new)

El cual tiene como centro unas 240 a 260 personas, lo cuál calza con la cantidad real indicada, ya que esta se encuentra entre estos valores y en el centro de la distribución.

Tercera parte, probabilidad de mordeduras, dada una mordedura

# Dataset nuevo
df.mordidos = df[df$bites_month_1 == 1,]
# Printeamos la cantidad de carteros mordidos el primer mes, y también
# la cantidad de carteros mordidos el segundo mes dado que fueron mordidos
print(sum(df.mordidos$bites_month_2))
## [1] 56
print(sum(df$bites_month_1))
## [1] 146
new <- rbinom(1e4, size=sum(df$bites_month_1), prob=samples)
simplehist(new, xlab= "new bites predictions")

Claramente este valor se encuentra fuera del centro teórico mostrado por el gráfico de densidad del predictive. En este caso, el modelo no es útil para realizar consultas indexadas ni tampoco considerar subconjuntos del dataset, debido a que la creación del posterior y de los samples considera a ambos meses de forma conjunta.

Pregunta 3: Inferencia Sobre Dos Parámetros

En esta pregunta se trabajara con el dataset “notas.csv” el cual contiene las notas históricas de un curso desconocido. Suponga que los datos vienen de una distribución \(\mathcal{N}(\mu,\sigma^2)\), el objetivo de la pregunta es estudiar el comportamiento de los datos y los posibles valores de \(\mu,\sigma\) mediante técnicas de inferencia Bayesiana.

Usted sabe un dato extra sobre la información, los valores de \(\sigma\) en la grilla se mueven en el intervalo \([0.5,1.5]\) y, además, tiene una fuerte creencia a que es mas probable encontrar la desviación estándar real entre \([0.5,1]\) que en \((1,1.5]\). De hecho, estudios señalan que la probabilidad de encontrar sigma en los valores \([0.5,1]\) es de 2/3, mientras que 1/3 para el resto de intervalos.

  • Modifique el siguiente código que permite realizar una grid approximation para \(2\) parámetros. Proponga priors para \(\mu\) y \(\sigma\), justifique su elección.
# Leer información
data_notas <- read.csv("notas.csv")

# Función para crear likelihood dado mu y sigma
grid_function <- function(mu,sigma){
   prod(dnorm(data_notas$Notas, mean = mu, sd = sigma))
}

# Valores de la grilla
grid_mu <- seq(from=1, to=7, length.out=100)
grid_sigma <- seq(from=0.5, to=1.5, length.out=100)

# Se crea la grilla 2d
data_grid <- expand_grid(grid_mu,grid_sigma)

# Se guarda la likelihod
data_grid$likelihood <- map2(data_grid$grid_mu,data_grid$grid_sigma, grid_function)

# Se transforma el forma de map2 a una columna
data_grid <- unnest(data_grid,cols = c("likelihood"))

# Valores de los priors

# Como no tenemos información a priori de la media, asumiremos que todos son equiprobables
prior_mu <- rep(1 , 100) 
# Por otro lado sabemos que la primera mitad de valores es el doble de probable que la otra mitad
prior_sigma <-  c(rep(2/3 , 50), rep(1/3, 50))

# Se crea la grilla 2d de priors
prior <- expand_grid(prior_mu,prior_sigma)

# Se calculan los valores del prior
data_grid$prior <-  map2(prior$prior_mu,prior$prior_sigma, prod)
data_grid <- unnest(data_grid,cols = c("prior"))

# Se calcula el posterior
data_grid$unstd_posterior <-  data_grid$likelihood * data_grid$prior

# Se estandariza el posterior
data_grid$posterior <- data_grid$unstd_posterior/sum(data_grid$unstd_posterior)

# Se ajustan los valores de la posterior para que no sean valores tan pequñeos
data_grid$posterior <- (data_grid$posterior - min(data_grid$posterior))/(max(data_grid$posterior)-min(data_grid$posterior))
  • Tras haber ejecutado el código de la parte anterior ejecute el siguiente, ¿Que puede decir de los valores de la distribución? Se recomienda hacer modificaciones en el código para realizar un mejor análisis y estudio.
# Punto de referencia
# Se recomienda cambiar estos valores por unos adecuados que le permitan estudiar
# Los valores de la distribución de mejor manera

row <- data_grid[which.max(data_grid$posterior),]
valor_x <- row$grid_mu
valor_y <- row$grid_sigma

# Grafico

punto_comparacion <- tibble(x = valor_x, y = valor_y)

plt <- data_grid %>%
  ggplot(aes(x = grid_mu, y = grid_sigma)) +
  geom_raster(aes(fill = posterior),
    interpolate = T
  )+ 
  geom_point(x = valor_x, y = valor_y, size = 1.3,color="white")+
  geom_label(
    data = punto_comparacion, aes(x, y),
    label = "Punto Comparación",
    fill = "green",
    color = "black",
    nudge_y = 0, # Este parametro desplaza la caja por el eje y
    nudge_x = 1 # Este parametro desplaza la caja por el eje x
  )+
  scale_fill_viridis_c() +
  labs(
    title = "Posterior para Mean y Standard Deviation",
    x = expression(mu ["Mean"]),
    y = expression(sigma ["Standar Deviation"])
  ) +
  theme(panel.grid = element_blank())

plt

Con la distribución anterior podemos decir el el valor más probable se encuentra con la media entre 6 y 6.2, mientras que la desviación estandar entre 0,5 y 0,75. Es importante recalcar que esta distribución está en tres dimensiones, por lo que es posible dibujar una superficie en 3d donde se muestre la curva de distribución con el peak en el punto de comparación.

  • A continuación se presenta un código que permite realizara la distribución dada por sampling from grid approximation ¿Para que sirve este proceso? ¿Que puede deducir del gráfico?
# Codificamos los datos
x <- 1:length(data_grid$posterior)

# Sampleamos los indices
posterior_samples_aux <- sample(x,size = 1e4, replace = T, prob = data_grid$posterior)

# Obtenemos los verdaderos valores de la sampling distribution
posterior_samples <- data_grid[posterior_samples_aux,]

# Obtenemos solos los valores relevantes para la densidad
df <- data.frame(posterior_samples$grid_mu,posterior_samples$grid_sigma)

# Realizamos las densidades
dens(df)

Lo anterior es útil para poder visualizar de forma separada la distribución de la posterior por cada uno de los ejes, que en este caso es la media y el otro la desviación estándar. Con el gráfico podemos ver claramente que los valores de mayor probabilidad de la media están entre los valores 5.9 y 6.2, mientras que para los valroes de la desviación estándar, se encuentra entre 0.5 y 0.7.

Es claro que para distintos valores de prior estas distribuciones cambiarían, en especial para la desviación estándar, en donde no asumimos equiprobabilidad entre valores de los rangos, si no que la primera mitad es el doble de probable que la segunda mitad de la distribución.

 

A work by CC6104